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ABSTRACT 

With appropriate modifications, the Hoffman-Ribak algorithm that constructs con- 
strained realizations of Gaussian random fields having the correct ensemble properties 
can also be used to construct constrained realizations of those non-Gaussian random 
fields that are obtained by transformations of an underlying Gaussian field. For example, 
constrained realizations of lognormal, generalized Rayleigh, and chi-squared fields hav- 
ing n degrees of freedom constructed this way will have the correct ensemble properties. 
The lognormal field is considered in detail. 

For reconstructing Gaussian random fields, constrained realization techniques are 
similar to reconstructions obtained using minimum variance techniques. A comparison 
of this constrained realization approach with minimum variance, Wiener filter recon- 
struction techniques, in the context of lognormal random fields, is also included. The 
resulting prescriptions for constructing constrained realizations as well as minimum vari- 
ance reconstructions of lognormal random fields are useful for reconstructing masked 
regions in galaxy catalogues on smaller scales than previously possible, for assessing the 
statistical significance of small-scale features in the microwave background radiation, 
and for generating certain non-Gaussian initial conditions for TV-body simulations. 

Key words: methods: data analysis - large-scale structure of Universe. 



INTRODUCTION 



C/3 



The recovery of a signal from noisy and incomplete data is 



a classic inversion problem tliat is of particular relevance 
when analysing observations of the galaxy distribution for 
two reasons. First, regions like the 'zone of avoidance' that 
is a result of obscuration by the Galactic Plane prohibit the 
construction of complete all-sky catalogues. The existence of 
'masked' regions raises intriguing questions about the con- 
nectivity of structures separated by a masked region, and so 
about the scale of the largest coherent structures in the ob- 
served Universe. Secondly, if one assumes that the luminous 
galaxies that are observed sample an underlying density field 
that is smooth, then the discreteness of objects introduces 
'shot-noise' that may inhibit the reconstruction of the un- 
derlying smooth field from the data. 

Recently, Lahav et al. (1994) addressed these problems 
of incomplete sky coverage and shot-noise within the con- 
text of whole-sky galaxy surveys. They note that, when the 
underlying density field is a Gaussian random field, mini- 
mum variance techniques (e.g. Rybicki fe Press 1992) and 
the conditional probability framework of constrained realiza- 
tions (Hoffman fe Ribak 1991) both provide the same esti- 
mate for the 'optimal reconstruction' of the underlying field, 
given the observed field. This is essentially a consequence of 



three properties peculiar to Gaussian random fields. First, 
a Gaussian field is completely specified by two parameters, 
essentially its mean and variance, so that minimizing the 
variance between an observed Gaussian field and the re- 
constructed Gaussian field will certainly produce an opti- 
mal estimate of the underlying field. Secondly, for Gaussian 
fields, the most probable field is also the mean field, so that 
reconstruction of the mean field, given the observed field, 
also gives the optimal estimate of the true underlying field. 
Thirdly, the distribution of residuals around the mean real- 
ization of a Gaussian random field is also Gaussian, with a 
variance that is independent of the value of the mean real- 
ization. As Hoffman fe Ribak (1991) showed, with the pre- 
vious two properties, this third property admits a powerful, 
efficient technique for constructing constrained realizations 
(e.g. Bertschinger 1987) of Gaussian random fields. 

Counts-in-cells analyses of the angular distribution of 
IRAS galaxies show that, on large angular scales (angu- 
lar separations 9^ 10°), the /i?j4 5 distribution is essentially 
Gaussian (e.g. Sheth, Mo fe Saslaw 1994). Therefore, when 
Lahav et al. (1994) apply their minimum variance technique 
and construct constrained realization reconstructions of the 
angular distribution of IRAS galaxies, to recover structure 
on scales ^15°, their results are spectacular. On smaller 
scales, however, the IRAS distribution is extremely non- 
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Gaussian. Unfortunately, when the underlying distribution 
is non-Gaussian, the mean field reconstruction and the most 
probable reconstruction are not, generally, the same. It may 
be that reconstructions that minimize variance will differ 
significantly from constrained realization techniques in this 
regime. 

Sheth et al. (f994) show that, on nonlinear scales (cor- 
responding to small angular separations), the data appear 
consistent with the hypothesis that the IRAS galaxies are 
drawn from an underlying smooth density field that is a 
lognormal. Bernardeau fe Kofman (f995) argue that there 
is good reason to expect a lognormal to provide a good ap- 
proximation to the true density field traced by galaxies, if 
the initial power spectrum was similar to that which is ex- 
pected in the standard cold dark matter model. Therefore, 
to extend the Lahav et al. (f 994) analysis to smaller, nonlin- 
ear scales, a lognormal model for the non-Gaussian density 
field is likely to be both useful and plausible. This is be- 
cause, for a lognormal distribution, it is possible to compute 
the difference between the minimum variance estimator of 
the underlying field, the mean field, and the most probable 
field, given the observed field. A systematic study of the dif- 
ferences between these estimators of the true underlying log- 
normal field will allow optimal whole-sky reconstructions of 
the galaxy density distributions even on the very small scales 
where the galaxy distribution is highly nonlinear. Moreover, 
it will help quantify the strengths, and also the limitations, 
of the minimum variance and the constrained realization ap- 
proaches to reconstructing an underlying density field. 

In the cosmological context, the ability to construct con- 
strained realizations of random fields is also useful because 
it allows the generation of initial conditions for numerical 
simulations that are designed to incorporate, a priori, par- 
ticular features of interest. For instance, rare density peaks 
of a Gaussian random field may behave quite differently from 
those in non-Gaussian fields. The study of the detailed dy- 
namics of the collapse of density peaks is greatly facilitated 
by the ability to 'make them to order' in the initial con- 
ditions of a numerical simulation. Moreover, the problem 
of 'cosmic variance' due to the long-wave modes that are 
missing in a finite-sized simulation box may be partially al- 
leviated by using those long-wave modes measured observa- 
tionally to constrain the initial conditions in the simulation. 

Section 2.1 summarizes useful background material 
about Gaussian and lognormal random fields. The notation 
follows that of Coles fe Jones (f99f). Section 2.2 computes 
the conditional lognormal distribution; it is a 'shifted' log- 
normal. Then it computes the distribution function of the 
residuals of a lognormal field, given a set of constraints that, 
say, specify the value of the field in certain regions (so the 
distribution of the constraints is also lognormal). Since, for 
lognormal fields, the mean and most-probable fields are dif- 
ferent, the distribution of residuals from both fields is cal- 
culated. For both cases, the distribution of the residual field 
depends on the values of the constraints. In this respect, 
lognormal fields are different from Gaussian fields. 

It would seem, therefore, that the constrained realiza- 
tion algorithm suggested by Hoffman fe Ribak (f99f) can- 
not be applied to lognormal density fields. Finally, Section 
2.3 demonstrates that the Hoffman-Ribak algorithm can be 
used to construct constrained realizations of the log of a log- 
normal density field. It also shows that taking the exponen- 



tial of the constrained realization of the log of the lognormal 
field gives a density field that is the 'optimal' constrained 
realization. That is, when used in this way, the Hoffman- 
Ribak algorithm can construct an ensemble of constrained 
realizations of lognormal fields that has the correct proper- 
ties, just as it can for Gaussian fields. Section 3 provides 
an algorithm that produces constrained realizations of log- 
normal fields that can be compared directly with, e.g., the 
IRAS catalogue, while Section 4 relates all these results to 
the use of minimum variance, Wiener Filter reconstructions 
of lognormal random fields. This last is of interest given the 
demonstration (Rybicki fe Press f992; Lahav et al. f994) of 
the equivalence of the constrained realization and minimum 
variance approaches when reconstructing Gaussian random 
fields. 

Section 5 discusses extensions of the Hoffman-Ribak al- 
gorithm to treat other random fields that are obtained by 
simple transformations of an underlying Gaussian field. It 
argues that, analogous to the lognormal case, an appropri- 
ately modified application of the Hoffman-Ribak algorithm 
will also produce constrained realizations of these other non- 
Gaussian random fields, and that these constrained realiza- 
tions will possess the correct ensemble properties. As spe- 
cific examples, it shows this to be true for chi-squared fields 
having ra-degrees of freedom, and for generalizations of the 
Rayleigh and Maxwell distributions considered by Coles fe 
Barrow (f987). 



2 CONSTRAINED REALIZATIONS OF 
LOGNORMAL RANDOM FIELDS 

2.1 Notation 

The one-point probability distribution function (pdf), fi{x) 
for a Gaussian random field, X{r), is 



fi{x)dx 



dx 



a/2^ 



exp 



{x - l^ f 

2(72 



(1) 



where fi and a are, respectively, the mean and the variance 
of X. For a Gaussian random field all the higher order n- 
point pdf's, f„{x), of field values at different positions r;, 
are multivariate Gaussians so that 



fn{x) = (27r)--|Mr-exp 



where x = (xi, . . . , x„ 
matrix with elements 



.(2) 



M, 



x{ri), and M is the covariance 



(3) 



where all the /ti's are assumed to have the same value, fi. 

The lognormal field (LN) is obtained simply by trans- 
forming a Gaussian field via 



Y(r) = exp [X(r)], 

so that its one-point pdf is 



fi{y) dy 



dy 1 
y a/2j 



exp 



2(72 



(4) 



(5) 
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where /t and cr^ are the mean and variance of the underlying 
Gaussian field X. The corresponding ra-point pdf is 



fniVl, 



(27r)-"''"|Mr'''" X 



TT -!- X exp 
r"--"- yi 



(6) 



where M is the covariance matrix of the X values defined 
above. 



2.2 Properties of the constrained field 

This subsection uses the notation of the simple example of 
Section 2 in the Lahav et al. (f 994) Wiener filter letter. The 
simple case they considered demonstrates all the interesting 
properties without having to introduce more complicated 
notation for a random field rather than a random variable. 

First, define the following statements: 

P(a;) = the probability the random variable has value x, 

P{x,y) = the joint probability that one variable has 
value X and the other is y, 

P{x\y) = the probability that one variable has the value 
X given that the other is y. 

To construct constrained realizations, think of x as the 
value of a realization given that the constraint has value y. 
Next, calculate the mean value of an ensemble of realizations 
of x, in which each realization of x is subject to the (same) 
constraint y, {x\y), and then calculate the distribution of the 
residual fluctuations around this value of the mean realiza- 
tion, i.e., the distribution of u = x — {x\y). The key to the 
Hoffman-Ribak algorithm is that, for a Gaussian random 
field, the distribution of u is independent of the value of y. 
Is this true for a lognormal? 

For clarity (and so that factors of fi don't proliferate) 
the following considers the variable x which is assumed to 
be distributed lognormally and for which the mean, fi, of 
the underlying Gaussian variable is zero. Then, equation (5) 
shows that 



f f 



exp 



(lux) 



2(ji 



(7) 



where cr^ = {(Inx)^). When /t = 0, equation (6) shows that 
the joint probability 



-PLN(a;,2/) 



IMI-- f 
■ exp 

l-Kxy 2 



MZ^{\ax) 



+2Mr2 (In x\ny)+ M^"/ (In y) 



phere M is the covariance matrix: 

.2 



M 



with a. 



(8) 



(9) 



{{\nxf), 



{(In 2/)^), and { = {lux Inj/). So, 



{ involves the cross-correlations between In x and In y, the 



determinant |M 



is just (7^(7^ 



- { , and the inverse, M 



easily computed. Then, the conditional probability is 



1 1 



X exp 



2 \ aWy - e 



In X 



{In 2/ 



X A/27i 



exp 



(In X 



(10) 



where the final expression defines /t' and cr'^. Recall that 
when X and y are each Gaussian distributed, the conditional 
distribution of a; given j/ is a shifted Gaussian. Equation (fO) 
shows the analogous result for a lognormal distribution: 
when X and y are each distributed lognormally, the con- 
ditional distribution is a lognormal distribution (compare 
equations fO and 5), whose associated underlying Gaussian 
is shifted. (The mean and variance of the underlying shifted 
Gaussian are /t' and cr'^, respectively.) This means that if 
X and y are distributed lognormally, then an ensemble of 
realizations of x, with each realization subject to the (same) 
constraint y, will have a lognormal distribution. 

Since the mean and the most probable values of a log- 
normal distribution are different, this means that the most 
probable value of x given the constraint y, does not equal 
{x\y), the mean value of x given y. This illustrates one 
important difference between lognormal and Gaussian dis- 
tributions. When constructing constrained realizations of a 
lognormal field, it may be important to decide whether to 
construct realizations that resemble the mean realization, or 
whether to reconstruct the most probable realization. Below, 
expressions for both possibilities are derived. 

The most probable value of x given a value of y is that 
value of X for which 



(^Pi.n(x\y) 
Ax 

so that 



0, 



{x\y) 



exp 



ilny 



(11) 



where /t' and cr'^ were defined in equation (fO). On the other 
hand, the mean value 



{x\v)i_ 



d In a; e'' 
\/2i 



X exp 



<y'i<y% 

2 



In X 



ilny 



exp 



ilny 



2<J% 



.m' + (<t'"/2) 



(12) 



When { ^ (i.e., the limit where lux and In j/ are not 
correlated), then (a;|2/)mp a^mp, and {x\y)^^ {x)i,uj 
as one would expect. In the other extreme, when x and 



y are completely correlated, so that { 



then 
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(a;|2/)mp = {^\y)i,]si = 2/) which is also expected. It is impor- 
tant to notice that, in general, equations (11) and (12) may 
be quite different from each other. 

Next, it is necessary to compute the distribution func- 
tion of the residuals from the mean realization (and also 
from the most-probable realization). To compute the dis- 
tribution function of the residuals from the mean define 
n = X — {x\y)j^j^ and compute P(m|2/). Then find the mean 
and the variance of this distribution. Recall that when x 
and y are Gaussian random variables, the mean and second 
moments of u are independent of y. Here, 



(13) 



(18) 



(n) = J nP(n\y)dn = 
as expected. However, 



'P{u\y)du = {x"|2/> - {x\yf = e"'''+"'^" - 6""' + "" , (14) 



so that, as equations (10-12) show, the distribution of the 
residual field depends on the value of y. 

Similarly, define the residuals from the most probable 
realization v = x — (a;|2/)mp and compute the distribution 
function P{v\y). Then 

{v\y) = {x\y) - {x\y)rap, (15) 
and 

{v^\y) = {x^\y) - 2 {x\y) {x\y)rap + ix\y)l,p, (16) 
so that the distribution of v also depends on y. 



2.3 Generalization to random fields 

Although the analysis above only considered the random 
variables x and y for which the mean of the underlying Gaus- 
sian variable was zero, it is clear that when considering a 
random field rather than a random variable, with some non- 
zero value for the underlying mean Gaussian field, nothing 
significant will be different. For this generalization to log- 
normal random fields the mean constrained realization is 



{/(^)|r>. 



exp 



^6(r')«,^Mn< 



(17) 



where the lognormal field is f{r) = /(fi, . . . , r„), T is the 
set of, say, m constraints each of which specify, for example, 
that the value of the (lognormal) field at the ith point is c;, 
£,i{r) is the cross correlation between the underlying Gaus- 
sian field (that is associated with f{r), the lognormal field) 
at r and the log of the ith constraint, {(0) is the variance 
of the underlying Gaussian field, and fij is the matrix that 
specifies the correlation between the log of the ith and j'th 
constraints (it is the covariance matrix of the underlying 
Gaussian field that is associated with the lognormal con- 
straint field). (This notation is similar to that of Hoffman fe 
Ribak 1991.) 

Similarly, the most probable realization is 



(/(^)|r) 



exp 



and the variance of the residuals around the mean realiza- 
tion, F{r) = f{r) — {/(?") |r)LN, is easily calculated from 
equation (17). As with the simpler random variable exam- 
ple considered earlier, this generalization to a random field 
shows that when the density field and the constraint field are 
both lognormal, then the distribution functions of the resid- 
uals from both the mean and the most probable realizations 
are a function of the constraints. 

The results above suggest that, since the variance of 
the residuals is a function of the constraints, the Hoffman- 
Ribak technique may not be used to construct constrained 
realizations of lognormal fields. However, because a lognor- 
mal field is so easily related to its underlying Gaussian field 
(equation 4), there is a simple transformation that simplifies 
the problem considerably. 

Equations (5) and (6) suggest defining a new E field 
so that e = Iny. Then the one-point pdf of E is Gaussian 
and, as equation (6) shows, all the formalism developed for 
Gaussian random fields may be applied to the E field. For 
example, the ra-point pdfs of the E field are all multivariate 
Gaussians. Now, equation (10) shows that the conditional 
distribution of a lognormal field, when the constraints are 
distributed lognormally, is just a lognormal whose underly- 
ing Gaussian is shifted. Therefore, for the E field defined 
here, the conditional distribution for the E field is just this 
shifted Gaussian (since, if the constraints are lognormal, 
then the log of the constraints is Gaussian). Therefore, by 
judicious choice of the constraint field, the Hoffman-Ribak 
algorithm may be used to construct constrained realizations 
oiE. 

In particular, if the constraint field essentially specifies 
the value of e (that is, it specifies the value of the log of 
the density field y) at a given set of points, then, if is a 
lognormal, (so E is Gaussian), the distribution of the con- 
straints on Y will be lognormal, so that the distribution of 
the constraints on E will be Gaussian. Since the E field and 
the constraint field are both Gaussian, and the conditional 
distribution is just a shifted Gaussian, the distribution of 
the residuals (from the mean of E) in any particular real- 
ization of the constrained field will be independent of the 
values of the constraints. So, the Hoffman-Ribak algorithm 
can be applied to construct constrained realizations of E. 
Then, simply taking the exponential of E at every point 
gives a constrained realization of the lognormal field. The 
question is: how optimal is this procedure for constructing 
constrained realizations of the lognormal field, Y? 

The following example illustrates this question. Assume 
that a constrained realization of the (Gaussian) E field 
has been constructed (using the Hoffman-Ribak algorithm). 
This means that in between the points where the field is con- 
strained to have certain values the field fluctuates around 
its mean realization value, with most of the fluctuations oc- 
curring within 'one standard deviation' or so. Since E is 
Gaussian, this is perfectly acceptable as a realization of E. 
However, we are really interested in = exp E. Taking 
the exponential of the E field, at the points where the con- 
straints have been specified, gives a Y field that (by con- 
struction) has the correct values at these points. In between 
these points, however, as a result of the exponential trans- 
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formation, the 'one standard deviation' departures from the 
mean E realization become much larger departures from the 
mean Y realization. It is not obvious that these fluctuations 
correctly sample the true (lognormal, cf. equation 10) dis- 
tribution of fluctuations around the mean (lognormal) Y 
realization. This section shows explicitly that this proce- 
dure (i.e., transforming the constrained realization of the 
underlying Gaussian field, E) does, indeed, give the opti- 
mal procedure for constructing constrained realizations of 
lognormal fields. 

The mean and variance of the underlying constrained 
Gaussian variable E are fi' = ^\ny/ay and cr'^ = (cr^cr^ — 
{^)/(7^, respectively. The mean and the most-probable val- 
ues of the associated lognormal variable Y are '^^'^ 
and e** ""^ . Comparing with equations (11) and (12), we 
see that (a;|2/)mp and {a;|2/)LN have precisely the same values 
as the mean and most probable values of the Y field that 
is obtained by taking the exponential of the constrained E 
field. 

Furthermore, let U (R) denote the variance of the resid- 
uals from the mean Y = exp E realization. That is, define 



U{R) 



(19) 



The first term in the angle brackets is Y, written explic- 
itly as the exponential of the value of the variable E, which 
is written as the sum of its mean value, /t', plus a resid- 
ual, R, from that mean. The second term is the mean 
value of the lognormal variable Y. Now, E is Gaussian dis- 
tributed, with mean fi' and variance cr'^. This means that 
the residuals, R, from the mean E realization have a Gaus- 
sian distribution with mean zero, and variance cr'^, so that 
p(R) = ''^'^ /\/2Tra'^ . This, with equation (19), implies 
that 



U{R) 



2ii' + 2R 



■ 2e 



M' + fl„M' + (<T'=/2) 



2e 



+ e 

<t'2/2 



e^p(R) dR 



' J P(R) dR I 



2/Ll' 



2/Li +2ct ■ 



2e" 

2fl -t-cr 



72 <t'^/2 , a 

' e + e 



= e — e 

where we have used the fact that 



(20) 



nR 



e"^p{R) dR 



nR e 



-R^/2<t'- 



V27r(7'2 



■d_R 



-( R-ncj'^ f /2o 



\f2-h 



■dR : 



(21) 



Comparing equations (20) and (14) shows that the vari- 
ance of the residuals around the exponential of a constrained 
Gaussian variable, U{R), is the same as the variance of the 
residuals around the mean of a constrained lognormal vari- 
able (equation 14). 

Thus, this calculation shows that the ensemble of log- 
normal Y fields that are obtained by taking exponentials 
of an ensemble of (appropriately) constrained realizations 
of E fields, correctly samples the underlying constrained 



lognormal distribution of realizations. In other words, in- 
sofar as the Hoffman-Ribak algorithm provides an optimal 
technique for constructing constrained realizations of Gaus- 
sian random fields, the technique here, of constructing con- 
strained realizations of the log of a lognormal field, and then 
taking the exponential of the constrained realization, is op- 
timal for constructing constrained realizations of lognormal 
random fields. Thus, although the direct calculations of the 
previous subsection (and the complications involved in equa- 
tions 14-16) suggested that the Hoffman-Ribak algorithm 
would not be optimal for constructing constrained realiza- 
tions of a lognormal random field, with the transformation 
described here, the Hoffman-Ribak algorithm can, indeed, 
be used to construct optimal constrained realizations of a 
given lognormal density field. 

To illustrate the method. Fig. 1 shows the steps in- 
volved in constructing a constrained realization of a log- 
normal distribution, given one constraint, say Ci , which 
specifies the value of the field at the origin. The top panel 
shows a cut through an unconstrained realization of a three- 
dimensional Gaussian random field that has power spectrum 
P{k} oc exp[— (fci?^)^]. The Gaussian field is generated 
on a 32^ grid using FFT techniques, and is similar to the 
one presented in the top panel of Fig. 1 in Hoffman fe Ribak 
(1991). The mean realization, given the value at the ori- 
gin, /(O), is shown by the dotted line. It is determined by 
/(r) = /(0){(r)/{(0), where {(r) = {(|r - 0|) is the auto- 
correlation function (i.e., {(r) is the Fourier transform of 
the power spectrum), and {(0) is the variance of the field. 
The dashed lines show the mean realization plus and mi- 
nus one standard deviation, given the value at the origin, 
where the standard deviation (7(r) = -\/{(0) - K(»-)V?(0)] 
(cf. Hoffman fe Ribak 1991). At every point r, the residual 
is obtained by subtracting the value of the mean realization 
at r (dotted line) from the actual value of the unconstrained 
field (solid line) at that position. 

The dotted line in the middle panel shows the mean 
realization of the Gaussian when it is constrained to have a 
2{(0) peak at the origin, (i.e. /(O) = Ci = 2{(0)): f(r) = 
Ci {(f)/{(0), and the dashed lines show the mean realization 
plus and minus one standard deviation. The solid line shows 
this mean plus the residual calculated from the previous 
panel; it is the constrained realization of the 2{(0) peak. 

The solid line in the bottom panel shows the constrained 
realization of the corresponding lognormal field; it is ob- 
tained simply by taking the exponential of the solid line in 
the second panel. The dotted line shows the exponential of 
the mean of the constrained Gaussian realization. The dot- 
dashed line shows the ensemble mean of the constrained 
lognormal realizations, and at each point r, it is calculated 
using /(f) = exp[/t'-|- ((7'^/2)], where fi' = Ci {(f)/{(0), and 
cr'^ = {(0) — [{('')^/{(0)]. The triple dot-dashed line shows 
the most-probable values of the constrained lognormal field, 
and is given by /(f) = exp[/t' — cr'^]. These relations gener- 
alize those derived in equations (11-16). 

The remainder of this section contains a brief digression 
on the philosophy of the approach of dealing with the log of a 
variable rather than its actual value. The next section gives 
a prescription for constructing a constrained realization of, 
for example, the IRAS catalogue. 

Although the analysis above deals mainly with build- 
ing constrained realizations of density fields, it can easily be 
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Figure 1. Constrained realizations of Gaussian and lognormal 
fields. (Top) Solid line shows a cut through an unconstrained 
3-dimensional Gaussian random field that has power spectrum 
P{k) oc exp[— (fci?^)^]. Dotted line shows the mean value 

of the realization given the value at the origin; dashed lines 
show mean plus and minus one standard deviation. (Middle) 
Solid curve shows the constrained Gaussian field constructed 
from the previous panel: it is constrained to have a peak that 
is two-standard deviations high at the origin. Dotted and dashed 
lines show the mean constrained realization, and the mean plus 
and minus one standard deviation, respectively. (Bottom) Solid 
curve shows the corresponding constrained lognormal field. Dot- 
ted curve shows the exponential of the mean field of the underly- 
ing Gaussian, i.e. the exponential of the dotted curve of the pre- 
vious panel. Dot-dashed curve shows the mean constrained log- 
normal realization, and triple dot-dashed curve shows the most- 
probable constrained lognormal realization. 



generalized. For example, it may be that one person (called 
A) is interested in building constrained realizations of some 
density field as measured by, say, its optical flux, whereas 
someone else (called B) is more interested in the field as 
measured by its optical magnitude, rather than by its op- 
tical flux. The latter is essentially the log of the former. 
If B were to build a constrained realization of the magni- 
tude distribution, A could take the appropriate exponential 
transformation of B's realization to obtain an estimate of the 
distribution of flux. Provided the distribution of magnitudes 
is Gaussian, so that the distribution of flux is lognormal, not 
only will taking the exponential of the constrained realiza- 
tion of the magnitude distribution give a good indication of 
what the optical flux distribution looks like, but, as shown 



above, A's exponential transformation of B's constrained re- 
alization provides an optimal constrained realization of the 
distribution of flux. 

Returning to the problem of constrained realizations of 
density fields, this suggests that if the true density distri- 
bution is really lognormal, then the convenient variables in 
which to work are not the values of the density field at a 
given point, but rather, the log of those values. If one is in- 
terested in just the log of the density field, and the density 
field is a lognormal, then the HofFman-Ribak algorithm will 
give the 'optimal' constrained realization of that lognormal 
field. Furthermore, comparison of this constrained realiza- 
tion with the log of, say, the IRAS distribution, is a fair way 
of comparing the constrained realization to the data. Per- 
haps the way to put this is as follows: If the density field, 
as measured by its IRAS flux is lognormal, then, when ap- 
plied to the log of the flux, the HofFman-Ribak algorithm 
will give the optimal constrained realization of the magni- 
tude distribution. Taking the exponential of the constrained 
realization will give an optimal realization of the flux field. 



3 CONSTRUCTING THE CONSTRAINED 
FIELD 

The best way to construct a constrained realization of a 
lognormal field that is similar to, say, the IRAS distribu- 
tion, even on small highly nonlinear scales, is as follows. 
Construct a constrained realization of the E field using the 
HofFman-Ribak algorithm. Next, obtain the Y field by tak- 
ing the exponential of the E field. Then use the prescription 
of Coles and Jones (1991, Section 8.3) to convert the density 
field Y into a point distribution. One could even arrange the 
conversion so that the number of points in the constrained 
realization is the same as the number of galaxies in the IRAS 
catalogue. Finally, smooth the point distribution with a filter 
of choice to look at the density distribution on any scale of 
interest. Compare with the result when the IRAS catalogue 
itself is smoothed with that same filter. 

The all important first step is to construct a constrained 
realization of the E field. This requires knowledge of the co- 
variance function of the log of the IRAS distribution; once 
this covariance function is known, the following steps are 
straightforward. This covariance function is obtained as fol- 
lows. The two-point correlation function of the random vari- 
able y is 



{{yi - {y))iy2 - (y))) 
{y? 



(22) 



where r = \ri — ^2], 2/i and 2/2 are the values of the field at 
ri and ^2, and {«/) is the mean value of y. If y is distributed 
lognormally, then its mean can be calculated similarly to 
equation (12); it is 



{y) = exp \ n + 



(23) 



where fj, and cr^ denote the mean and variance of the under- 
lying Gaussian. For this lognormal, {(r) is easily calculated 
from the two-point pdf: 



!+«('•) 



1 

w 



2/12/2 /(2/1, 2/2)d2/i Ay2 
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1 tI ' 



(24) 



where (t\ = Mn = M22 = <y\, and p = M12 = ^(In yi — 
2/2 — /(*))• This final term shows that the covariance 
function, p, is related to the correlation function of the un- 
derlying Gaussian (the correlation function is p/ p? ). Writing 
H(r) for the covariance function of the Gaussian means that 
p = H(r), so equations (23) and (24) show that 



l + {(f) = exp [H(f)]. 



(25) 



(e.g. Vanmarcke 1983; Coles fe Jones 1991). 

As Coles fe Jones note, this exact relation resembles the 
Politzer fe Wise (1984) approximation used in biased clus- 
tering theory. Politzer fe Wise were trying to calculate the 
correlation function for high-level regions of a Gaussian ran- 
dom field. They found that, to a good approximation, the 
correlation function, {(f), of high-level regions of a Gaus- 
sian field is related to the correlation function, H(r), of the 
underlying Gaussian by 1 -|-{(r) = exp[H(r)]. Since galaxies 
are thought to be biased traces of the matter distribution, 
this provides another (albeit heuristic) justification for con- 
sidering a lognormal field as a good model for the galaxy 
distribution — the correlation function of the high density 
regions of a Gaussian field (the dark matter distribution) 
looks just like the correlation function of a lognormal field 
(the observed galaxy distribution). 

Equation (25) shows that, although it has not been mea- 
sured directly, the covariance function of the log of the IRAS 
distribution, H(r), can be obtained from the IRAS correla- 
tion function itself. So, to construct a constrained realiza- 
tion of the E field only requires a reliable measurement of 
the /i?j4 5 correlation function, {(f), and oi p = {In y). 

To make constrained realizations of the projected dis- 
tribution, {(r) should be replaced by io(d). The IRAS an- 
gular correlation function can be obtained in a number of 
ways: (1) directly measuring io(d), (2) via Limber's equa- 
tion from the spatial correlation function measured by, e.g., 
Saunders, Rowan-Robinson fe Lawrence (1992), (3) from the 
3-D power spectrum (Fisher et al. 1993), or (4) from the an- 
gular power spectrum. Although Fourier transforming the 
IRAS power spectrum is the natural choice when making 
constrained realizations of the large scale distribution (e.g. 
Lahav et al. 1994), since this paper is concerned with smaller 
scale features (hence the move from the Gaussian to the log- 
normal), it is more natural to use the correlation function 
itself, since, on these small scales, it is a better estimator 
of the correlations. Most measurements of the correlation 
function [using either (1) or (2) above] give 



(2G) 



with 7 K, 1.6 and ^0 ~ 0.11° out to about 6° after which it 
drops quickly to zero. 

With equation (26) for the correlation function, and 
with the relations given earlier that relate the covariance 
function of the underlying Gaussian to the correlation func- 
tion of the lognormal, it should be relatively straightforward 
to construct constrained realizations of lognormal random 
fields that are similar to the IRAS distribution. 



4 MINIMUM VARIANCE 

RECONSTRUCTIONS OF LOGNORMAL 
RANDOM FIELDS 

If a quantity whose true value is s is measured to have the 
value y = s -\- n, where n is some noise on the measure- 
ment, then the Wiener filtered minimum variance linear re- 
construction of s, given the measured y, is 

5wf = S[S-FN]-' y = Fy (27) 

where S = l^ss^), and N = {nn^). This relation holds what- 
ever the true distribution of s, whatever the measured dis- 
tribution of y, and whatever the distribution of the noise, n: 
there is no requirement that any of these fields be Gaussian 
(cf. Rybicki fe Press 1992). 

Lahav et al. (1994) show that if s and n are Gaussian 
random fields, so y is also a Gaussian random field, then 
the minimum variance reconstruction, Smv, is the same as 
the Wiener filtered reconstruction, which is the same as the 
most probable reconstruction defined using the constrained 
realizations formalism, and recall that, for Gaussian random 
fields, this most probable realization is also the mean real- 
ization. Since, in general, the mean and the most probable 
fields are different for a lognormal, it is interesting to com- 
pare the Wiener Filter reconstructed field, subject to the 
condition that the true field is lognormal, with the mean 
and most probable constrained lognormal fields. 

We consider two cases. First, we consider a lognor- 
mal signal s, in which the noise on the underlying Gaus- 
sian field is assumed to be Gaussian and additive, so that 
logy = log 5 -|- log Ti, where y is the measured value of s, 
and log s and log n are both Gaussian. (This means that 
y = sn; the noise is multiplicative, not additive.) Then it is 
straightforward to apply the minimum variance reconstruc- 
tion to this underlying Gaussian field (rather than to the 
lognormal field directly). Consider the field that is obtained 
by taking the exponential of the minimum variance recon- 
struction of the underlying Gaussian field: 



„(log «)wf 



exp F log y 



(28) 



Now, as mentioned above, the Wiener filter, minimum vari- 
ance estimator of a Gaussian field is also the mean of the 
underlying constrained field (e.g., Lahav et al. 1994). So, in 
the notation of the previous sections, this reconstructed field 
is simply , where p' denotes the mean value of the con- 
strained field. However, recall that the most probable value 
is e** , whereas the mean value is e** '^'^'^ (compare 
equations 11 and 12). Thus, the estimator given in equa- 
tion (28), i.e., the exponential of the Wiener filtered value 
of the underlying Gaussian field, has the attractive property 
of always being between the most probable and the mean 
constrained realizations. 

The second case we consider is when s is lognormal and 
n is Gaussian. This case of Gaussian noise, whatever the 
underlying signal, may be useful for reconstruction analyses 
of small-scale features in, e.g., the microwave background. A 
simple example that uses the random variable, rather than 
the random field will illustrate the results. Let y = s -\- n, 
where s is a lognormal distributed signal, and n is Gaussian 
noise. Then log j/ = log(s -\- n) = logs -|- log(l -|- n/s). In 
the limit where s is large relative to the noise, this simpli- 
fies to log y Si log s + (n/s). Since n is Gaussian distributed. 



8 Ravi K. Sheth 



for sufficiently large s, we can treat n/s as Gaussian dis- 
tributed also. Therefore, in this approximation, log y is the 
sum of two Gaussian variables, so it too is Gaussian. So, 
this case is similar to that discussed above. This means that 
exp (log s)wf = exp (_F log j/), is a good estimator of the un- 
derlying signal, since it is always between the most probable 
and the mean reconstructions. However, when the noise is 
not small, which is the more general and useful case, ob- 
taining a useful estimator of the underlying signal is more 
complicated. 

To obtain an expression for the most probable value of 
the signal given the measured value y, consider P{s\y) = 
P{s,y)/P{y) = P{s)P{y\s)/ P{y), where the distribution 
P{s) is lognormal, and P{y\s) = P{s + n\s) = P{n) is Gaus- 
sian. Since P{y) = J P{s)P{n}ds, with n = y — s, is merely 
a normalizing factor, there is no need to compute it explic- 
itly. Below, we will compare the minimum variance (Wiener 
filtered) estimate, Swf, of the underlying signal given the 
measured value y, with the most probable value of P{s\y). 
This most probable value, Smp, is given by requiring that 
dP{s\y)/ds = 0. This requirement means that 

dP{s\y} 1 d 



ds 



P{y) ds 



P(s)P(n=y-s) 



P{n=y-s) dP{s) P{s) dP{n=y-s) 
~~P{y) d^ + PW 

P{n}P{s} ( 1 lns\ P{s)P{n) {y - s) 



P{y) V s sal) P{y 
P{n)P{s) f y - s 1 Ins 



0, 



(29) 



where we have assumed that s is distributed lognormally 
in such a way that the mean of the underlying Gaussian, 
/is = (Ins) = and ^(Ins)^^ = (t1 , and that the noise 
is Gaussian distributed around zero with variance cr^. This 
implies that Smp and y satisfy the nonlinear relation 



V 



(j'i In s 



mp 



(30) 



On the other hand, the relation between the Weiner filter 
estimate, Swf, and the measured value, is linear. Namely, 



{{^-{^)f) 

{s - {s)f) + {n^) 



(31) 



Thus, the relation between Smp and s^t is complicated and, 
in general, the Wiener filter estimate differs from the most 
probable one: s^t ^ Smp. The difference between the two 
quantities depends on the variance, and hence on the power 
spectra of the signal and the noise. When the noise is vanish- 
ingly small, so that crj 0, then both Smp and Swf y, as 
they should. When the variance of the noise is much larger 
than that of the signal, then Smp e~'^' as expected (cf. 
equation 11). In contrast, the Wiener filter estimate is 0, in 
this limit. Finally, note that if s is distributed lognormally, 
then s must always be positive. Notice that the estimator 
Srap is positive even when y is negative. In contrast, s^t has 
the same sign as the measured value, y, so it can be negative. 

Since the relation between Smp and s^t is nonlinear, 
the accuracy of the Wiener filter estimate of the underlying 
signal is most easily illustrated numerically. Before doing so, 
it is useful to extend equation (30) to the case of a measured 



random field y = s -\- n, where s is an underlying lognormal 
signal and n is the noise associated with the measurement, 
and is assumed to be Gaussian. Then, the most probable 
estimator of the underlying signal s satisfies the relation 



5 -|- diag (^s S ^ In5-|-N '"5 = N ^ y, 



(32) 



where JV and H are the covariance matrices of the noise 
and of the log of the signal, respectively. Compared to the 
Wiener filter estimator (equation 27), this equation is much 
more difficult to solve. 

In the absence of correlations, both H and JV are diag- 
onal, so that the most probable estimator of 5 at a given 
point in the field is independent of the values of 5 or n at 
other points in the field. Thus, in the absence of correla- 
tions, equation (32) reduces to equation (30) at each point. 
Similarly, the Wiener filter estimate (equation 27) reduces to 
equation (31) at each point. Of course, this is consistent with 
equation (25) which shows that if the underlying Gaussian 
field is uncorrelated, so that its covariance matrix is diago- 
nal, then the covariance matrix of the associated lognormal 
field is also diagonal. 

Fig. 2 shows a numerical example of the difference be- 
tween the most probable and the Wiener Filter estimators 
of the underlying signal for the case when the noise, n is 
Gaussian and the signal, s, is lognormal, and both have di- 
agonal covariance matrices. The light solid line in the top 
panel shows y = s -\- n, whereas the heavy solid line shows 
the actual value of the signal, s. Both s and y are plotted in 
units of the root mean square value of the (Gaussian) noise. 
In the Figure, Srms = 5.23 at each point, which corresponds 
to an underlying Gaussian with zero mean and dispersion 
(Tg = 1.75. This value was chosen because it gives approxi- 
mately the same value of the variance of the galaxy density 
field on the scales on which the lognormal may be good ap- 
proximation. For convenience, the (Gaussian) noise has zero 
mean and dispersion equal to that of the signal, Spnis- 

Since the covariance matrices are diagonal, the Wiener 
filter estimator is s^t = yj'^ at every point. The light solid 
line in the bottom panel shows this Wiener Filter estimator. 
The heavy solid line in the bottom panel shows 5mp; at each 
point it is the largest root of equation (30) that is less than 
the measured value y at that point. For large values of y, 5mp 
is nearly the same as y, which means that it, in the Figure, is 
a factor of two larger than the Wiener Filter estimator (large 
values of y are all those that are greater than, say, 3(7n). For 
those values of y that are intermediate (say, between Cn 
and 3(7n), the most probable estimator is comparable to the 
Wiener filter estimator. For those values of y that are less 
than (7n, Smp is approximately e » ~ 0. Finally, note that 
whereas 5wf is often negative, 5mp is always positive. 

Fig. 2 shows that it is possible for 5wf to differ sub- 
stantially from 5mp. As stated above, the exact difference 
between Swf and Smp will depend on the details of the co- 
variance matrices of the noise and the signal (equations 27 
and 32). Thus, this Section shows explicitly that if the un- 
derlying signal is non-Gaussian, then the Wiener Filter esti- 
mator of the underlying signal, given a noisy measurement 
of that signal, may not be very close to the most probable 
value of that signal, given the same measurement. 
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Figure 2. Example of the difference between the most probable (heavy solid line, bottom panel) and the Wiener filter (light solid line, 
bottom panel) estimators of the underlying signal for the case when the noise, n is Gaussian and the signal, s, is lognormal, and both 
have diagonal covariance matrices. Light solid line in the top panel shows y = s -\- n. Heavy solid line in top panel shows the actual value 
of the signal, s. 



5 DISCUSSION AND EXTENSIONS 

It is straightforward to extend the technique developed in 
this paper for constructing constrained realizations of log- 
normal random fields to other non-Gaussian fields that can 
be obtained by transforming an underlying Gaussian field. 
Here we consider two other examples, the generalized Chi- 
squared random field with n degrees of freedom (hereafter 
Xn), and the generalized Rayleigh distribution; these ran- 
dom fields generalize some of the distributions considered 
by Coles fe Barrow (f987). 

The characteristic function of a p-variate x^-field is 



det||I-i2MT| 



-1/2 



(33) 



where G„{t) denotes G„{ti,t2, . . . ,tp), I is the identity 
matrix, T is the diagonal matrix with diagonal elements 
ti,t2, . . . ,tp, and M is the covariance matrix of the un- 
derlying Gaussian field (e.g. Lukacs f977). Although, in 
general, the corresponding probability density /n(z) = 



fn{zi, Z2, 



with 



(34) 



where the xjs are independent normal variates is compli- 
cated, a calculation of the case when p = 2 and the xjs have 
the same unit variance and normal distribution around mean 
zero, will serve to illustrate our purpose. Then, 

„2,. , . ^\ r(f) 



fn{zi,Z2) 
fn(zi) f„{z2) 



exp 




phere 



fu{z)Az : 



j("-2)/2 

T(fr 



Az 



(36) 



provided z > Q, and fn{z) = otherwise, T{z) is the Gamma 
function, Ia{z) is the modified Bessel function, and p = 
{X1X2) is the covariance function of the underlying Gaussian 
variables (e.g. Vere-Jones f967). 

From these expressions, and using the series expansion 
of the modified Bessel function, the conditional distribution, 
fn{z2\zi) = fn(zi,Z2)l f„(zi), is easily computed; it is 



fn(z2\zi)Az2 



^(n-2)/2g(z+A)/2 
2^ 



E 



\^z^ 



2dz2 

22'=fc!r(f -Ffc) (f-/92) 



(37) 



where z = 2z2/{l - f?), A = 2f?zil(l - f?). Equation (37) 
for the conditional distribution of Z2 given zi is a non- 
central Xn distribution (e.g. Miller f964, p. 56; Stuart fe 
Ord f99f. Sections 23.4, 23.6), so it can be understood as 
arising from a shifted Gaussian distribution. More generally, 
if /ti and /t2, and cr^^ and cr^^, are, respectively, the means 
and variances of the original Gaussians, the mean and vari- 
ance of the underlying (shifted) Gaussian distribution are 
n' = p-\/{zi/n}/al^, and cr'^ = {al^al^ - p'^)/al^, respec- 
tively. By analogy with the lognormal example discussed 
earlier (compare equation fO), it is clear that taking the ap- 
propriate transformation of the (appropriately) constrained 
Gaussian field will give constrained realizations of x^-fields 
that have the correct ensemble properties. 

Recently Bunn et al. (f994) used constrained realiza- 
tions of Gaussian random fields to assess the statistical sig- 
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nificance of hot and cold spots in maps of the microwave 
background radiation as seen by the COBE DMR. As Luo 
(1995) and Coulson et al. (1994) note, the data presently 
available admit the possibility that on small (~ 0.5°) an- 
gular scales, fluctuations in the radiation may have a non- 
Gaussian spatial distribution. Luo (1994) argues that Xn- 
fields provide a family of random fields that range from 
the highly non-Gaussian (low-ra) to the nearly Gaussian 
(large ra), and so they are able to model a large range of 
non-Gaussian distributions. In addition, temperature fluc- 
tuations from topological and non-topological defects in the 
framework of the 0(A^)(7-model are thought to be well de- 
scribed by a x^-field (Turok fe Spergel 1991). Both these 
considerations justify a continued interest in x^-fields. They 
also suggest that this generalization of the Hoffman-Ribak 
algorithm to x^-fields may be useful for assessing the statis- 
tical significance of structure in the microwave background 
on small angular scales. 

Additionally, as Sheth et al. (1994) show, the discrete 
Negative Binomial distribution provides a good fit to the 
distribution function of galaxy counts-in-cells. The Negative 
Binomial can be understood as arising from a Poisson sam- 
pling process [exactly the same Poisson process discussed 
by Layzer (1956) and Peebles (1980), and analogous to that 
used by Coles fe Jones (1991) to construct a discrete number 
count model from a (continuous) lognormal density field], 
applied to an underlying Gamma distribution which has the 
form shown in equation (36) (e.g. Stuart fe Ord 1991, Sec- 
tion 16). So, it appears that Poisson sampling of an appro- 
priately constrained realization of a x^-field, should produce 
constrained realizations of the Negative Binomial distribu- 
tion that are similar to the observed galaxy distribution. 

Coles Barrow (1987) studied the lognormal, the Xrn 
and two other non-Gaussian distributions that can be ob- 
tained by transforming a Gaussian, as models of the distri- 
bution of extreme hot spots in the microwave background. 
For completeness, we now study a generalization of their 
Rayleigh and Maxwell distributions. Expressions for the gen- 



eralized Rayleigh distribution, gn{i'i,r2, 



where 



\ 



(38) 



and where all the xjs are independent normal variates are 
given, e.g., in Miller (1964, ch. 2). Here, we simply consider 
the case when p = 2 [the case when p = 1 and ra = 2 is the 
Rayleigh distribution studied by Coles fe Barrow (1987), 
and they termed the p = 1 with ra = 3 case the Maxwell 
distribution]. Then, when the xjs all have zero means and 
unit variance. 



an(r) 
and 

an(ri,r2) 



2"/2r(f) ' 



(39) 



1 



(2p)(n-2)/2 (1 _ p2)r(._^ 



X exp 



0^] 

-{rf + rj) 
2(1 -p2) 



7(„ 



rir2P 
1 - p2 



(40) 



where p and T(r) are defined as for the Xn case. With these 
two expressions it is straightforward to show that g{r2\ri) 



is a Rayleigh distribution whose underlying Gaussian vari- 
ables are shifted [see Miller (1964) for the expression that 
describes a generalized Rayleigh distribution with non-zero 
mean and non-unit variance]. In general, the mean and vari- 
ance of the shifted Gaussians are fi' = pri / ^/na'^^ and 
cr'^ = {cr'^-^a'^^ — p'^j/cr'i-^, respectively. As shown explicitly 
for the lognormal case, it is clear that taking the appropri- 
ate transformation of the appropriately constrained Gaus- 
sian field will give constrained realizations of generalized 
Rayleigh fields that have the correct ensemble properties. 

All these results strongly support the speculation that, 
with the appropriate transformations of the constraint field, 
the Hoffman-Ribak algorithm is optimal for constructing 
ensembles of constrained realizations of all random fields 
that are obtained by transforming Gaussian fields. A sketch 
of the proof is given in the Appendix. 



6 CONCLUSIONS 

Section 2.3 showed that, with appropriate modifications, the 
Hoffman-Ribak algorithm can be used to construct an en- 
semble of constrained realizations of lognormal density fields 
that has the correct ensemble properties. Fig. 1 showed the 
steps involved in constructing such a constrained lognor- 
mal field. Section 3 provided an algorithm that can produce 
constrained realizations of lognormal fields that can be com- 
pared directly with, e.g., the spatial distribution of galaxies 
in the IRAS catalogue. 

Since the mean and the most probable values of any log- 
normal field are different, the accuracy of minimum variance 
reconstruction estimates of lognormal random fields is less 
than when reconstructing Gaussian random fields (for which 
the mean and the most probable are the same). Section 4 
provided two models for the relation between a measure- 
ment, the noise and an underlying lognormal signal. When 
the signal is lognormal and the noise is both lognormal and 
multiplicative (rather than additive), then the Wiener filter 
reconstruction provides a robust estimate of the underlying 
signal, since it always lies between the mean and the most 
probable values. However, when the signal, s, is lognormal 
and the noise, n, is Gaussian and additive, (so the measured 
quantity is y = s + n), then the accuracy of the Wiener fil- 
ter reconstruction depends on the power spectrum of the 
signal and of the noise. Equations (27) and (32) showed 
that, in general, the most probable reconstructed value is 
different from the Wiener filter estimate. Fig. 2 provided 
a simple numerical example of the way in which these two 
estimators differ. It showed that the most probable recon- 
struction differs significantly from the Wiener filter estimate 
at those points where the measured values are greater than 
3(7n, where is the root mean square of the noise. There- 
fore, when the underlying signal is lognormal, then the ac- 
curacy of the Wiener filter reconstruction will differ for dif- 
ferent data sets. 

Finally, Section 5 showed that the Hoffman-Ribak al- 
gorithm can be extended to construct ensembles having the 
correct ensemble properties of constrained realizations of all 
random fields that are obtained by transforming Gaussian 
fields. The proof was sketched in the Appendix. 
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APPENDIX A: CONDITIONED 
NON-GAUSSIAN RANDOM FIELDS 

Following, e.g.. Miller (1964), let X = {xi, X2, ■ ■ ■ , Xn} be 
an ra-dimensional Gaussian vector with the positive definite, 
symmetric, covariance matrix M. So 

f{X) dX = ^ ^ exp {-\x^M.-^X^ , (Al) 



(27r)"/2 



where f{X) dX denotes the probability that the vector lies 
between X and X + dX, and |M| denotes the determinant 
of M (compare equation 2). The covariance matrix, M, may 
be partitioned into apxp matrix, Ci, in the top left corner, 
an (n — p) X (n — p) matrix, C2, in the bottom right corner, 
a p X (n — p) matrix, B in the top right corner, and the 
transpose of B, B"^, an (n — p) x p matrix, in the bottom 
left corner. Since M is symmetric, Ci = C^ and C2 = C2 . 



Further, since M is nonsingular, its inverse, M"S exists, 
and it too can be partitioned. 

Let Q"'^, P"'^, R, and R"^ denote the corresponding 
partitions of M"'^. Since M is symmetric, Q = Q"^ and P = 
P^. Then PR^ = -B^C"^ and Q"^ = C"^ + RPR^ are 
Shur's identities, and |P| = |M||Ci|~'^, where the vertical 
bars denote the determinants of the matrices, is a special 
case of Jacobi's theorem (cf. Miller 1964, Section 1.3). 

Finally, let Xi = {xi, X2, ■ ■ ■ , Xp], with 1 < p < ra, and 
let X2 = {xp^i, Xp-\-2, . . . , Xn}- So, Xi is the p-dimensional 
vector comprised of the first p elements of X, and X2 is 
the (ra — p)-dimensional vector comprised of the last (ra — p) 
elements of X. This means that 

X'^M.'^X = XiQ~'^Xi + X2'R-'^Xi 

+X'[RX2 + x'^p-'X2 
= (X2 + PR^Xi)^p-'(X2 + PR^Xi) 

+ X?'(Q-' - RPR^)Xi, (A2) 

where the second expression follows by adding and subtract- 
ing Xi RPR Xi to complete the square in the X2 vari- 
ables. However, the second of Shur's identities shows that 



Xi(Q~^ - RPR^)Xi = XiC^^Xi. 



(A3) 



PR^ Xi means that 
-Vfp-\X2-V) + XlC^'Xi.{A4:) 
■PR^Xi = B'^C'^Xi, where the second 



So, setting V = — 

X^M~'^X = {X2 

Notice that V = - 
equality follows from the first of Shur's identities (compare 
with the terms in equation 17), is independent of X2- 



Now, it is clear that 
dXi 



/(^i)dXi 



exp 



— -^1 
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and recall that f{X) = f{Xi,X2)- So, the conditional dis 
tribution, f{X)/ f{Xi), is given by: 

f(Xi,X2)dXidX2 



f{X2\Xl)dX2 



f(Xi)dX^ 
dX2 



(27r)(' 
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Since |P| = |M||Ci|~'^ (Jacobi's theorem) and since V is 
independent of _X"2, this shows that ii f(Xi, X2) is Gaussian 
distributed, then the conditional distribution, f(X2\Xi), is 
a shifted Gaussian. 

Now consider the vector Y, where Y = g{X), and 
where J^f is a Gaussian vector with covariance matrix M. 
That is, the vector Y is some transformation of the Gaussian 
vector X. For example, if Y = exp X = {e^^ , e^^ , ■ ■ ■ , e^" }, 
then dX = dlogY = Yli dyi/yi, and the distribution 
of Y is lognormal (compare equation 6). Similarly, if 
Y = {xi, X2, ■ ■ ■ , x'^}, then y is a x^-(with one degree of 
freedom)-vector, and dX = d\/Y = Y[ <iyi/2^/yi (compare 
equation 37). Clearly, if Y = g{X), where X is a Gaussian 
vector, and ifY is partitioned into Yi and Y2, then the con- 
ditional distribution, f(Y2\Yi) = f(Yi,Y2)/f(Yi), will 
be given by the appropriate transformation of the underly- 
ing Gaussian that has been shifted hy V = — PR"^Xi = 
B^ C^^ (Yi ). As a result, with the appropriate modifi- 
cations described in the text, the Hoffman-Ribak algorithm 
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may be used to construct constrained realizations of Y that 
have the correct ensemble properties. 

Of course, this assumes that g~'^{Yi) is sinele valued. 
For the lognormal distribution this is satisfied. For the Chi- 
squared distribution, however, both the positive and the neg- 
ative square roots are valid solutions. In practise, this will 
not be a problem if the mean of the underlying Gaussian 
field is known to be positive, and if the (transformed) con- 
straints are unlikely to be negative. This will almost always 
be the case if, for example, the constraints are chosen to cor- 
respond to sufficiently high density peaks of the Chi-squared 
field. 



